%% analyze run-along-corridor task

% channels
%  4: forward speed
%  8: rotation speed
%  9: x-position (along corridor)
% 10: heading (2.5 is straight)
% 11: reward
% 12: y-position(across corridor)

animalID = 1
%animal ID 1:3 = M090 M091 M092

x_length = [10 30 60 100]
% x_length gives the length of the virtual corridor as an array, sorted by
% days

beh(1).M(1).day(1).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T9972.lvd');
beh(1).M(2).day(1).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T9973.lvd');
beh(1).M(3).day(1).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T9974.lvd');
beh(1).M(1).day(2).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10032.lvd');
beh(1).M(2).day(2).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10041.lvd');
beh(1).M(3).day(2).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10050.lvd');
beh(1).M(1).day(3).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10106.lvd');
beh(1).M(2).day(3).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10109.lvd');
beh(1).M(3).day(3).aux_data = load_lvd('D:\tempData\Behavior M090 - M092\S1-T10111.lvd');

for lnd = 1:length(beh)
    for knd = 1:length(beh(lnd).M)
        for gnd = 1:length(beh(lnd).M(knd).day)
            for ind = 1:length(beh(lnd).M(knd).day(gnd).aux_data(10,:))
                if beh(lnd).M(knd).day(gnd).aux_data(10,ind)>2.5
                   beh(lnd).M(knd).day(gnd).angle(ind)=(beh(lnd).M(knd).day(gnd).aux_data(10,ind)-2.5)/2.5*180;
                end
                if beh(lnd).M(knd).day(gnd).aux_data(10,ind)<2.5
                   beh(lnd).M(knd).day(gnd).angle(ind)=(abs(beh(lnd).M(knd).day(gnd).aux_data(10,ind)-2.5))/2.5*(-180);
                end
            end
        end
    end
end

figure,
id = 0;
for ld = 1:length(beh)
    for kd = 1:length(beh(ld).M)
        for md = 1:length(beh(ld).M(kd).day)
            id = id+1;
            subplot(length(beh(ld).M),length(beh(ld).M(kd).day),id), hist((beh(ld).M(kd).day(md).angle),100)
            xlim([-185 185])
            ylim([0 400000])
        end
    end
end

for lnd = 1:length(beh)
    for knd = 1:length(beh(lnd).M)
        for mnd = 1:length(beh(lnd).M(knd).day)
            beh(lnd).M(knd).day(mnd).polar = beh(lnd).M(knd).day(mnd).angle * pi/180;
        end
    end
end

for lnd = 1:length(beh)
    figure, rose(beh(lnd).M(3).day(1).polar,1000)
end

for ind = 1:length(beh)
    for knd = 1:length(beh(ind).M)
        for lnd = 1:length(beh(ind).M(knd).day)
            for mnd = 1:length(x_length)
                beh(ind).M(knd).day(lnd).x_length = x_length(1,lnd);
            end
        end
    end
end

fw_ch=4;
ro_ch=8;
xpos_ch=9;
ypos_ch=12;
ro_ch=10;
reward_ch=11;

for lnd = 1:length(beh)
    for knd = 1:length(beh(lnd).M)
        for gnd = 1:length(beh(lnd).M(knd).day)
            beh(lnd).M(knd).day(gnd).xpos = beh(lnd).M(knd).day(gnd).x_length*beh(lnd).M(knd).day(gnd).aux_data(xpos_ch,:);
            beh(lnd).M(knd).day(gnd).ypos = 2*beh(lnd).M(knd).day(gnd).aux_data(ypos_ch,:);
        end
    end
end

x_length = 25
aux_data = load_lvd('D:\tempData\S1-T13223.lvd');

xpos=x_length*aux_data(xpos_ch,:);
ypos=2*aux_data(ypos_ch,:);

vel_fw=[0 diff(aux_data(fw_ch,:))];
vel_fw(vel_fw>5)=vel_fw(vel_fw>5)-10;
vel_fw(vel_fw<-5)=vel_fw(vel_fw<-5)+10;
vel_fw=ftfil(vel_fw,1000,0,10);

vel_ro=[0 diff(aux_data(ro_ch,:))];
vel_ro(vel_ro>5)=vel_ro(vel_ro>5)-10;
vel_ro(vel_ro<-5)=vel_ro(vel_ro<-5)+10;
vel_ro=ftfil(vel_ro,1000,0,10);

for xnd = 1:length(aux_data(10,:))
    if aux_data(10,xnd)>2.5
       aux_data(10,xnd)=((aux_data(10,xnd)-2.5)/2.5*180);
    end
    if aux_data(10,xnd)<2.5
       aux_data(10,xnd)=((aux_data(10,xnd)-2.5)/2.5*180);
    end
end

fw_ch=4;
ro_ch=8;
xpos_ch=9;
ypos_ch=12;
ro_ch=10;
reward_ch=11;

aux_data = load_lvd('D:\tempData\S1-T13834.lvd');
reward_times=find(diff(round(aux_data(reward_ch,:)))==5);
nbr_rewards=length(reward_times)

heading_vec_x=smooth2(xpos-x_length*vel_fw.*cos(pi*(aux_data(ro_ch,:)-2.5)/2.5),10);
heading_vec_y=smooth2(ypos+x_length*vel_fw.*sin(pi*(aux_data(ro_ch,:)-2.5)/2.5),10);

figure;hist(aux_data(10,:),100)

for ind=60
    figure, hold on
    t_inds=reward_times(ind)+1500:reward_times(ind+1);
    plot(xpos(t_inds),ypos(t_inds),'linewidth',3);
    for knd=1:100:length(t_inds)
        plot([xpos(t_inds(knd)) heading_vec_x(t_inds(knd))],[ypos(t_inds(knd)) heading_vec_y(t_inds(knd))],'r');
    end
end

figure, subplot(3,2,6)
hold on
for ind=34
    t_inds=reward_times(ind)+1500:reward_times(ind+1);
    plot(xpos(t_inds),ypos(t_inds),'linewidth',3);
    for knd=1:100:length(t_inds)
        plot([xpos(t_inds(knd)) heading_vec_x(t_inds(knd))],[ypos(t_inds(knd)) heading_vec_y(t_inds(knd))],'r');
    end
end

